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ABSTRACT 


A study was made of the Space Shuttle Rocket Motor Casing during water 
impact. The problem was assumed to be static and equivalent static loads were 
used. The objective was to ascertain that current finite element analysis 
techniques could be relied on for design purposes. 

Preliminary analysis showed that 'the doubly curved triangular shell 
elements were too stiff for these shell structures. The doubly curved quadri- 
lateral shell elements were found to give much improved results. 

A total of six load cases were analyzed in this study. The load 
cases were either those resulting from a static test using reaction straps to 
simulate the drop conditions or under assumed hydrodynamic conditions re- 
sulting from a drop test. The latter hydrodynamic conditions were obtained 
through an emperical fit of available data. 

Results obtained from a linear analysis were found to be consistent 
with results obtained elsewhere with NASTRAN and BOSOR. 

The nonlinear analysis showed that the originally assumed loads would 
result in failure of the shell structures. The nonlinear analysis also showed 
that it was useful to apply internal pressure as a stabilizing influence on 
collapse. 

A final analysis with an updated estimate of load conditions resulted 
in linear behavior up to full load. 
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INTRODUCTION 
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In this report we describe the analysis of the Space Shuttle Rocket 
Motor Casing. The analysis was performed by the finite element method and 
use was made of doubly curved shell elements. The MARC-CDC program was used 
in this analysis. This program is described in the Appendix. 

The object of the analysis was to verify the use of such analysis 
by comparing its results with experiment. A further point of interest was 
to determine if the structure would behave linearly up to full load. The 
shell structure is loaded, in practice, .by inertia loads. A technique was 
developed during the study to account for the pseudo-static loads and 
reactions due to inertia. 

The analysis is divided into four parts: 

(a) In the initial stage the structure was modelled by means of 
triangular shell elements. These were found to give too stiff 
a result. (Sections 2 and 3) 

(b) Subsequently, a convergence study was made to compare the 
triangular shell element with the quadrilateral shell element. 
The quadrilateral element was found to be a much better element 
for this type of geometry and loading (Section 4). 

(c) The quadrilateral element was then used for analysis of the 
static model and the drop test model in the linear and nonlinear 
regime. (Sections 5 and 6) 

(d) Finally, a new estimate of the actual load during a drop test 
was used for an incremental analysis. (Section 7) 
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■ STATIC TEST MODEL ANALYSIS WITH CURVED TRIANGULAR SHELL ELEMENT 
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Idealization 

The static test model consists of a cylinder with two spherical closures 
at the ends. The model is shown in Figure 2.1. In the finite element idealization 
238 shell elements were used with a total of 160 nodes (Fig. 2.2). Appropriate 
boundary conditions were imposed to account for symmetry and remove rigid body 
motion. Analysis assumptions are discussed in the following sections. 


Thickness 

The assumed uniform thickness for the cylinder and the spheres are 0.347" 
and 0.231“ respectively. These thicknesses are based on the mean values employed 
in constructing actual experimental models. The stiffening effect at the joints 
is taken into account only at the end areas of the cylinder as shown in Figure 2.1. 
In these areas the shell elements are uniformly thickened consistent with the 
cross-sectional area of the joints. Uniform thickness of 0.444" and 0.534 are 
used corresponding to fwd and aft end joint areas respectively. 

It is noted that a beam-shell model was initially generated to simulate 
the stiffening effect. However, this model was abondoned because of the cost 
involved in the analysis of that model. 


Boundary Conditions 

Since the structure and loading are symmetrical with respect to the xz 
plane (Figure 2.1), the following boundary conditions are imposed. For all nodes 
on the xz plane: 




0 , 



0 
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In addition to the symmetry conditions the following boundary con- 
ditions are imposed to prevent rigid body motion: 


At Node A: 


u - 0 and w = 0 


At Node B: 


u = 0 

Note that these latter conditions should impose no loads. 


Loading v 

The assumed loading for the present analysis is shown in Figure 
2.3, and is defined as follows in terms of the surface coordinates e 1 (=Rx) 
and 6^ (=z). The pressure is varied linearly along the half length of the 
cylinder. Following the notation defined in Figure 2.3, the variation of peak 
section pressure P along the length (at =0) is: 


P = 


P.-P a 
b a 

T"T 

6 l” 6 _ 

b a 


(« - 


e a> + P a 


0 ) 


The pressure is varied also along the circumference as a cosine function 

so that the pressure vanishes at e^ t where 0 ^ is defined as a linear function 
2 

of 0 ^: 


e 1 = 




9 1 


0? - e' 




s a> +e 1 


( 2 ) 
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Therefore, the pressure p at a general point (e 1 , 8^) is given by 
following expression: 


«2\ _ r b ^a r „2\ , „ ~\ r7r i e 1 


P (e » e ) 55 [ -f-f (0 1 -e£) +P a ] . cos [Jj- ( TpTr 


■)] (3) 


b a 


44 <e 2 -4 + ®a 

e b' e a 


The numerical values assumed for the present analysis are: 


P a = 45. psi 

a 


P K = 8. psi 


I Peak values of pressure on symmetry plane 


“a = 15 ° < 5 1 = % • 

5 b = 60 ° ( 5 b ■ “b • R > 


circumferential extents of pressure 
loading 


The experimental strap support system was assumed to be frictionless, 
so that it provides a uniform normal pressure over the entire strap angle 
(120° - 60° in the symmetric half model). Then this reactive pressure, p^ as seen 
in Figure 2.4, was computed as a function of axial position (e 2 ) to provide 
section equilibrium in the x-direction, on the assumption that this was a 
reasonable modelling of the experimental procedure. 

Section equilibrium is achieved as follows: 

The component of the applied load, F^, at a section defined by e 2 is 

A 

given by the following: 

a a 

F® = / P . COSa . R . da (4) 
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where p Is given by Eq. (3) and a is defined by * a . R . 
Rewriting Eq. (3) in terms of a and substituting into Eq. (4) one 

obtains 


= R . P / COS a COS j “ da 


The above integration yields 


2ttR 
0 


F° = R . P . n o 

x (4) 2 - 4 

0 


. cos 


e 1 


(5) 


The component of the reactive load at a corresponding section is 


<f> 

F r = J p . cos 4 > . R d<j> = p . R sinj 
X o r 

Equilibrium condition Fj" = F^ gives the reactive pressure: 

A A 

2ttR 




P = 
r 


cos 


e 1 


sin| [(I R } 2 _ 4] 

e 


( 6 ) 


l is the strap angle used in the actual experiments and assumed to be 60°. 

In the MARC program the distributed loads are computed by numerical ' 
integration over each element. A subroutine defines the load magnitude point- 
by-point — seven points are used per element in the element used in the present 
analysis. To check the accuracy of this procedure, independent integrations 
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of the total applied loads in the x-direction (Equation (5) integrated with 
respect to 0 ) and in the y-direction were performed. These values were 
checked against the total distributed load in the x and y directions which were 
internally calculated by the computer program for use in the stress analysis. 
The comparison of the values is as follows: 

Numerical Integration 
Program 

Deviation (percent) 

RESULTS 

A small displacement elastic analysis was performed. The numerical 
results were obtained in terms of displacement plots as well as contours of 
direct and equivalent stresses on the top and bottom surfaces of the shell. 

Figure 2.5 shows the radial displacements on the symmetry plane along 
the length of the cylinder. In this analysis a maximum radial displacement of 
3.39" occurred at the mid-section. 

The cross-sectional deformation pattern is plotted in Figure 2.6 at 
the section where maximum radial displacement is obtained. 

Figure 2.7 shows the y component of displacement on the y-z plane. 

In the cylinder a maximum compressive hoop stress of 30,000 psi was 

t 

obtained on the surface close to where the maximum pressure was applied. The 
stress was mainly caused by bending. At the same point the compressive axial 




229,311 lbs. 
234,734 lbs. 


54,045 lbs 
54,391 lbs 


0.6 
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stress was 23,000 psi. It is noted, however, that a stress concentration 
occurs in the aft end sphere due to ovalization. A maximum compressive 
stress of about 51,000 psi was found in the hoop direction in the area (c) 
indicated in Figure 2.1. 




THICKNESS: 

CYLINDER 0.347" 

SPHERE 0.231" 

THICK SHELL ELM. FWD: 0.444" 

AF T •' 0.5 34" 


PROPERTY : 

E = 29 X 10 psi 
P-0 . 32 


BOUNDARY CONDITIONS 
SYMMETRY(XZ PLANE) 

v=0,-^- = 0 1 -^-*0,-i!S. 

se 2 se' $e' 


NODE A: 
u = 0,w*0, 

I 

Be 

NODE B: 


u = 0 


/FIG. e. I STATIC TEST MODEL 










•2744) (749.0744) ( 874 . 8744 ) 


(1000.6744) 

(1009.0544) (1066.6644) 
(1025.7644) 



U | (1161.8723) 

I ' (1130.562 1 ) 

' FIG. 2.2 (CONTINUED) 0 098.6644> 




L= 835.08 


FIG. 2.3 APPLIED PRESSURE 







$ =60° STRAP ANGLE 


FIG.2.4 REACTIVE PRESSURE 





FIG. 2:6 CROSS-SECTIONAL DEFORMATION AT MID-SECTION 





DROP TEST OR INERTIA LOADED MODEL 
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The geometric model was identical to that used for the strap load 
cases. This allows direct comparison between inertia and strap loads. 

LOADING 

The inertia loading is derived below.. This derivation assumes that 
the casing behaves essentially as a rigid body during the deceleration, that 
is, that for the computation of d'Alembert forces the distortion may be neg- 
lected. Then referring to Figure 3.1, let the accelerations of the casing at 
(x, y, z) be a x , a y , a z . From symmetry considerations and the rigid body 
assumption it follows immediately that 


at all points. 

The rigid body assumption gives additionally 

a x = A x + 2 e y 

and a - A - x e 

c y 


0 ) 


( 2 ) 


where A x , A 2 are the translational accelerations of the 
point (0,0,0) in the x- and z-directions , and is the 

rotary acceleration of the same point about the y-axis. 

* 

The three equilibrium equations remaining after symmetry con- 
siderations (about the x-z plane) are 
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(3) 


/ a pt dS = 0 
S z 

in the z-direction 

/ a pt dS = / p dS 
S 1 x 

3 s' 

in the x-di recti on 

and / (a z-a x) pt dS = / zp dS 

r A l i X 

S 1 

about the y-axis 


( 4 ) 


( 5 ) 


where s is the surface of the casing, 

is the surface subject to pressure, 

P x is the component of pressure in the x-direction 

1 2 1 o 

t (e , e ) is the shell thickness at a point (e , e ; 
and p is the shell density (assumed constant) 

Equations (2) and (3) give immediately 

A z = 0 (6) 

since t = t (z) is not a function of x. 

Equations (2) and (4) give 

pA / tdA + p 8 / zt dS = / p dS (7) 

x S y S $ 1 x 

’ and (2) and (5) give 

pA / zt dS + pe‘ /(z 2 +x 2 ) tdS = / zp dS (8) 

S y 1 x 

s 1 
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Equations (7) and (8) may be solved to give A x and P 6 y as 



(I y " 2c V 2 

tv- i y - V) 


( 9 ) 





2 

-V. I } 

y 


V) 


( 10 ) 


where P 

A 



dS = net external 


force in the x-direction 


z 

c 


/ z P v dS = 

d X 


center of influence of external pressure 


V = / tdS = volume of shell 


and 


J = / zt dS = 1st moment of area of shell about the 
y S y-axis 

I = / (z 2 +x 2 ) tdS = 2nd moment of area of shell about 
y S the y-axis. 


The load case is then defined by 

■jo 1 

p (e , 0 ) on S (as used in the previous analyses) 

and d'Alembert forces -pta in the x-direction, -pta in the 
z-direction per unit surface area at all points of the shell. 

The forces -pta x , -pta 2 are available by combination of equations 
(9) and (10) with (2). 
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In practice the constants defined after equations (10) were computed 
numerically by suitable coding in the shell element in MARC, resulting in the 
following values: 

P v = -2.347 x 10 5 lb 

z P = -1.6207 x 10 8 lb-in 

C X 

V = 7.20881 x 10 4 in 3 

0 = 3.569 x 10 7 in 4 

1 = 2.44714 x 10 10 in 5 

Then the coding was further modified to include -pta and -pta z 
as forces on the shell in the x- and z-directions. This technique (using nu- 
merical development of the inertia loads) made the inertia load analysis quite 
straightforward. The values of P , z , and V were checked by approximate hand- 

A L 

calculation to eliminate gross programming errors. 

RESULTS 

The results of this analysis are summarized in Figures 2-9. Figure 3.2 
shows the cross-sectional distortion at the most severely distorted section of 
the main cylinder. Figure 3.3 shows radial displacements on the x-z plane (the 
symmetry plane). In both cases the strap loaded case (static test case) results 
are shown for comparison. It may be seen that the inertia case always gives 
significantly less deformation. Comparison with the static model results 
shows rather similar stress distributions, but with reduced magnitude in the 
inertia load case. The equivalent stresses in the last ring of elements of the 
aft closure (the most highly stressed elements) are compared in Table 3.1 with the 
linear analysis results for the strap loaded case. 


3-4 



DISCUSSION 


In general the comparison with the strap loaded case indicates that, 
while the static test rig should predict the deformation mode with reasonable 
accuracy, the magnitudes of the displacements in such a test are 30% too high 
for the actual inertia loads and peak stresses are 60% too high. This is in- 
tuitively explained by noting that the net reaction from the strap includes a 
significant y-direction component not present in the inertia load. This should 
tend to increase the stress levels, but reduce the distortion. More significantly, 
the inertia load is generated uniformly around the circumference, while the strap 
load is confined to the lower of 120° segment. Thus, some of the pressure load is 
reacted directly by inertia loads, while in the strap loaded case all of the 
pressure load must be carried around to the opposite side of the shell. It is > 
this difference which seems most likely to give the large differences in stress 
and displacement magnitudes. The fact that the deformation modes are so similar 
suggests that the geometry renders the structure relatively insensitive to the 
distribution of particular load components, since in one case (straps) the load 
is confined to one section of the casing, while in the other case (inertia), 
while the net loads must be the same, the distribution is rather more uniform. 
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EQUIVALENT STRESS (psi) 




INSIDE 

OUTSIDE 


Element 

Strap 

Inertia 

Strap 

Inertia 

Number 

Reacted 

Reacted . 

Reacted 

Reacted 

225 

21592 

13543 

\ 

50711 

31525 

226 

18004 

12131 

16886 

10757 

227 

17552 

11261 

27211 

17407 

228 

39549 

24317 

12128 

7830 

229 

31012 

18718 

21071 

13722 

230 

15867 

9520 

22520 

13227 

231 

15635 

■ 9447 

15182 

9812 

232 

11227 

6977 

46572 

27973 

233 

32745 

20112 

41425 

24621 

234 

16288 

9411 

20792 

12106 

235 

17056 

10085 

24579 

14569 

236 

39022 

23514 

27150 

16558 

237 

18928 

11573 

36178 

22135 

238 

i 

14995 

8772 

18574 

11381 


TABLE 3.1 

EQUIVLANET STRESS 

IN END RING OF 



ELEMENTS ON AFT CLOSURE 







F1G.3.3 RADIAL DISPLACEMENT ON SYMMETRY PLANE 


MESH CONVERGENCE PROBLEM FOR SOLID ROCKET MOTOR 
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At this point in time the first results from the tests of the static 
model and parallel analysis at NASA with NASTRAN and BOSOR suggested that the 
results from MARC were four times too stiff. The NASTRAN analysis with 2000 
elements in an octant model predicted a mid-section displacement of 12.8 ins. 
in the static test and BOSOR predicted 13 ins. A mesh convergence study was 
then performed with a simplified model with the same governing parameters as 
the solid rocket motor booster to try to explain this discrepancy. 

In order to study the element efficiency for the actual rocket motor 
problem, a smaller, simpler analysis problem was investigated. This problem was 
felt to have the essential characteristics of the actual problem, but is of 
considerably smaller size, so that mesh convergence experiments were not prohi- 
bitively expensive. This section describes this smaller problem, and the results 
obtained. 


MODEL 

The model is one-eighth of a circular cylinder, with symmetry conditions 
on three edges and diaphragm support on the third (Figure 4.1). The loading is 
a normal pressure with the same distribution as that of the real problem (Figure 
4 . 2 ). 

1 2 

The pressure p at a general point (e , e ) is given by the following 
expression: 


P <e\ e 2 ) - [T~£(9 2 -9 2 ) +p,] 
e b e a a a 


cos [|- (— ^ 


0 


->] 


e b" e a / 2 2 V , -1 
"2 2 (e '"a* + 8 a 
e b' e a 
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The numerical values assumed for the present analysis are: 


P 

a 


45. psi 1 


peak values of pressure on symmetry plane 


= 8. psi 


-1 

a 

a 


-1 

% 


15° (e^ = 0g . R) 


60“ (e’ = ej . R) 


I circumferential extents of pressure loading 


1 2 

Where 9 and e are the Gaussian surface coordinates with 


1 2 
9 = Rot, 0 = z. 


Small displacement, linear elastic behavior was assumed. 


ELEMENTS 


Two thin-shell elements were used in this convergence test. Both use 
Koiter-Saunders shell theory, and both have C continuity (piecewise C ). Rigid 
body motions are represented exactly because of the isoparametric formulation 
of the elements. Loading was generated by consistant distribution. The tri- 
angular element is the Dupuis Element [1] which uses cubic interpolation 
corrected by rational functions. The quadrilateral element is a simple bi- 
cubic element using the nodal degrees of freedom 


i au 1 3U 1 


3 2 u^ 


1 2 * 12 

36 39 36 36 * 


i=l ,2,3 


at four nodes. This element is described by VerHague and Mallett in [2]. It 

should be noted that the element is only admissible as a rectangle in the 
1 2 

(e , 9 ) plane (the mapped plane of surface coordinates). This is not as severe 
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a restriction as might first be thought, since any surface coordinate set 
(e 1 , e 2 ) is allowable. However, the element is probably only useful for 
problems of the present type, where a regular mesh is satisfactory both in the 
Q^e 2 plane and on the real surface. For such problems, the element is ideal, 
since it is easy to form because of its straightforward interpolation functions, 
and is cheap to run, since nine (Gaussian) integration points in the (9 -9 ) 
plane prove sufficient. 

RESULTS 

In all cases regular grid spacings were used. The Dupuis element was 
used with five meshes: 4x6 divisions, 4x8, 6x8, 6x10 and 8x10, where the 

divisions around the quarter circumference are given first. The bicubic ele- 
ment was used with three meshes; 2x2, 3x3 and 4x4 elements. Further meshes 
were not run because it was felt little additional improvement would be ob- 
tained. The displacement results are characterized by the peak x-direction 
displacement (at the top center of the cylinder) in Figure 4.3. These indicate 
that any of the bicubic meshes gives an essentially converged solution, while 
the Dupuis triangle shows rather slow convergence. This came as a surprise, 
since the Dupuis element has shown rapid convergence in the classical cylindrical 
roof problem of Scordelis and Lo [3]. However, the relatively poor performance 
of the triangle was observed in some problems by Dupuis [3]. 

CONCLUSIONS 

The above mesh convergence study was taken to indicate that the Dupuis 
element was unsuitable for the rocket motor casing analysis, but that the bi- 
cubic element would provide reliable results with any reasonable mesh spacing. 
The convergence studies included small displacement, linear elastic behavior 
only, but because of the completeness of the element formulation (and its iso- 
parametric form), it was anticipated that it should give comparably good results 
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in the nonlinear range. Indeed, an axi symmetric form of this element was shown 
to have excellent convergence properties by McNamara [5]. 
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A BI-CUBIC QUADRILATERAL ELEMENT 
X DUPUIS TRIANGULAR ELEMENT 
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STATIC AND INERTIA LOAD COMPARISON - LINEAR ANALYSIS 
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With the results of the mesh convergence study in hand, it was de- 
cided to reanalyze the static and inertia loading cases. 

This section covers a linear comparison between static (strap- 
reacted) and inertia loadings of the 120" rocket motor casing. This analysis 
is preliminary to a similar comparison, based on non-linear incremental analysis. 
In all cases the pressure pad is on the aft half cylinder, with its peak at 
the mid-section of the cylinder. This linear comparison has been made before 
in section 3. The difference is that the present analysis uses a higher order 
quadrilateral shell element which has been shown to give rapidly convergent 
displacement solutions for this type of loading in section 4. 

MODEL 

The modelling is of the same 120" casing described in previous sections. 

The use of the new quadrilateral shell element allowed modelling with a coarser 
mesh at small penalty in solution accuracy, with considerable savings in computer 
cost. The geometric model assumed for the present analysis is shown in Figure 
5.1a, and the mesh used is shown in Figure 5.1b. 

Several points must be emphasized with regard to this geometric model- 
ling: 

1. Uniform thickness has been assumed throughout the cylinder. 

This is an unrealistic, but conservative assumption. The 
difficulty with the inclusion of thickening effects is that 
the thickened lengths are quite short compared to the element 
size used in the analysis. Because of the cost of mesh re- 
finement to capture thickening effects it was decided to use 
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the un-thickened model. Thickening could be achieved 
through the use of beam elements as stiffeners at the 
joints. While this approach is adequate for small dis- 
placement elastic analysis, stress predictions are not 
very good, so that this modelling is unsuitable for non- 
linear work. It must be emphasized that the un-thickened 
model will give conservative results, because the peak 
stresses occur at the aft (loaded) cylinder/sphere joint, 
which should be thickened. 

2. The aft spherical closure has been modelled as a complete 
hemisphere. Previous analyses by MARC have assumed a 60" 
opening in this closure, with no additional stiffness. Our 
analysis with the new quadrilateral shell element showed this 
end condition to be a most critical parameter, for both the 
peak displacements on the cylinder and the stressing of the 
closure itself. A strap-reacted analysis with the open ended, 
unstiffened closure showed severe ovalization of the sphere with 
associated high stresses. The closed end model was adopted 
after discussion with NASA personnel, since the actual casing 
would have some additional stiffness at this point from the 
nozzle and other members. The complete closure assumed here 
may possibly be over-stiff, in which case stresses on the clo- 
sure would be under-estimated. 

3. The mesh convergence tests indicated that a course 

mesh would be adequate for displacement predictions with 
the new shell element. Since the non-linear analysis 
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will be quite lengthy, it was decided to remesh with 5 
sections around the half model (180°) and 12 sections along 
the cylinder. The results reported in section 4 indicate 
that this mesh spacing will give reliable displacement pre- 
dictions for pressure pad loading. The linear analyses re- 
ported here took about 12 minutes each on a CDC 6600 - a 
considerable savings over previous analyses. 

LOADING 

The loadings used in the present analysis follow the same patterns 
reported previously. In both cases the pressure distribution is the same as that 
shown previously in Figure 2.4. Since the load is distributed by numerical in- 
tegration in the program, check sums are compared to an 'exact' integration. These 
sums give the net force in the x (vertical) and y (horizontal) direction due to 
the pressure on the half cylinder: their values are as follows: 


F x (lb) F(lb) 

•Exact' 229311 54045 

Previous Triangular 

element mesh 234734 54391 

Present Quadrilateral 

element mesh 228500 53793 


For the inertia case, rigid body acceleration fields, a x (x, y, z) 
were computed by the technique defined in section 3. The five constants 
defining the acceleration fields are again .computed numerically; the following 
table compares their values: 
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Previous Triangular 
Element Mesh 

Current Quadrilateral 
Element Mesh 


Total Vertical Load, 

-2.347 x 10 5 lb 

-2.285 x 10 5 lb. 


Moment of x Component of 

pressure about 0 , z . P 
y c x 

-1.6207 x 10 8 lb-in 

-1.5898 x 10 8 lb-in 


Volume of Casing, V 

7.2088 x 10 4 in 3 

7.0629 x 10 4 in 3 


First Moment of Casing 

about 0 , Jy 

3.569 x 10 7 in 4 

3.5128 x 10 7 in 4 


Second Moment of Casing 

about 0 , I 

10 5 

2.4471 x 10 ,u in 0 

2.4033 x 10 10 in 4 


RESULTS 

a) Displacements 



* 


Displacement modes are quite similar in the two loading cases, 
with rather different peak displacement magnitudes. The dis- 
placed profiles of the top (pressure loaded) and bottom edge 
are shown in Figure 5.2 (scaled to 45 psi peak pressure), and 
the distorted shape of the most severely deformed radial section 
is given in Figure 5.3. Mote the asymmetry about the plane x = 0: 
both the shape and magnitude of the vertical displacements 
(Figure 5.3) are quite different for the two edges. This would 
suggest that analyses based on the assumption of symmetry about 
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the plane x = 0 would overpredict the displacement and 
stress response under pressure loading since the reactive 
loading (inertia or strap) produces a less distortion on 
the opposite side of the shell. That is, symmetry about 
x = 0 would be an overly conservative assumption. 

Peak displacement magnitudes, scaled to 45 psi peak pressure, 
are as follows: 

Strap Reacted Inertia Reacted 

Displacement at node 

61 17. 3" 13.6" 

An explanation of the reduction in displacement response 
was offered in [1]. This displacement is basically the same 
as the NASTRAN and BOSOR predictions. It should be noted that 
this analysis is for a half model as opposed to the quarter 
models used for the NASTRAN and BOSOR analysis, 

b) Stress 

Peak stresses now occur at the intersection of the cylinder 
with the aft spherical closure. This is in contrast to the 
previous comparison. [1], where peak stresses occurred on the 
spherical closure, because of the ovalization allowed by the 
unstiffened hole in the aft sphere of that analysis. Peak 
stresses, scaled to 45 psi peak pressure, predicted by the 
present analysis are as follows: 

Strap Reacted , Inertia Reacted 


Equivalent Stress (psi) 

87440 

68420 

Hoop Stress (psi) 

59530 

46860 

Axial Stress (psi) 

53400 

34550 
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NONLINEAR ANALYSIS FOR INERTIA LOAD WITH AND 

WITHOUT INTERNAL PRESSURE 6 


This section describes the results of two non-linear incremental 
analyses of the 120" rocket motor casing. The difference between the two 
analyses lies in the loading and geometry of the model. The loading pattern 
in the first analysis is the same as analyzed previously (section 3), but 
applied incrementally using initially steps of 10% of full load (40 psi peak 
pressure). The second analysis uses an additional uniform internal pressure 
of 2 psi together with incremental loading of the external load pattern, in 
5% steps. 

MODEL 

The geometric model assumed for the non-linear analysis with no internal 
pressure was identical to that analyzed previously described in section 3. For 
the analysis including internal pressure the model was modified so that the 
forward spherical cap was closed. The quadrilateral shell element described in 
section 4 was used in both models with the same mesh used in section 5. The 
model is shown in Figures 5.1 and 5.1a. 

It must be emphasized that the geometry used assumes constant thick- 
ness throughout the cylinder. The actual motor casing has thickened joints, so 
that the model is assumed to yield conservative results. 
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LOADING 


The loading patterns used in these non-linear analyses were identical 
to those used previously, with an inclusion of 2 psi internal pressure in the 
second analysis. Since the models are almost identical (the only difference 
is the inclusion of a closed forward sphere in the internal pressure case) the 
rigid body acceleration fields, a (x,y,z), remained essentially unchanged and 

A 

the values reported in section 3 were used for both analyses. 

RESULTS 

a) Displacements 

The displacements obtained from the two non-linear analyses 
are quite similar in mode because the external loading dominates 
the internal pressure. Figure 6.1 plots the distorted shape 
of the most severely deformed nodal section at maximum load 
reached. As was evident in the analysis of section 3 there is 
no symmetry about the plane x=0. 

Figure 6.2 and Figure 6.3 show the displaced profiles of the top 
(pressure loaded) and bottom edges at maximum load reached in the two 
load cases. Again the asymmetry is evident. Note, however, that although 
the radial displacements are very similar in both mode and magnitude, the 
results of the analysis run with 2 psi internal pressure cor- 
respond to 45% of full load, whereas the results of the analysis 
with no internal pressure are at 35% of full load. The solution 
does not converge to a specified tolerance after this load, in- 
dicating imminent collapse (stiffness singularity). Thus, the 
internal pressure is seen to stiffen the structure. This can be 
seen clearly on Figure 6.4 which presents load-displacement curves 
for node 61 (the most critical node). 
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b) Stress 

Peak stress magnitudes calculated by the MARC program occur 
in Element 41 and are as follows: 



Internal p=0 psi 

Internal p=2 


(35% of full load) 

(45% of full 

Equivalent Stress (psi) 

32,200 

50,000 

Hoop Stress (psi ) 

27,217 

52,300 

Axial Stress (psi) 

26,500 

36,611 


Stress components and the equivalent stress are contoured for both analyses in 
figures 6.5 - 16. These stresses are shown for the highest loads applied in each 
case. There is a high degree of similarity in the stress patterns in the two 
cases with the highest stresses occurring under the loading pad. 

Figure 6.17 plots stresses versus % load for the centroid of the highest 
stressed element, Element 41. Notice that for the 10% load incrementation case 
(no internal pressure) the solution in stress oscillates significantly between 
increments. For the 5% load incrementation analysis (2 psi internal pressure), 
the solution is much more stable until after 30% application of peak load. 

This is the usual behavior as the collapse load is reached. 

Recall that the internally pressurized case withstood 45% of full load 
as opposed to 35% in the non-pressurized analysis. This extra 10% of load ex- 
plains the sizeable differences in the peak stresses noted above. 

It must also be remembered that the results presented here are con- 
servative since there is no account taken of the thickening of the joints due 
to the spherical-cylinder intersection. 
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At this stage, new results from the drop test instrumentation re- 
sulted in a less severe pressure loading condition. The pressure was more 
localized. The analysis was repeated for this loading condition. 

A non-linear, inertia loaded analysis of the same 120" casing geometry 
analyzed in previous work has been performed. The pressure distribution dif- 
fered from that assumed previously, with a localized pressure spike of 55 psi 
at Z=450". Some comparisons with previous analyses using the same geometry 
but the previous pressure distribution are made in this section. 

MODEL 

The modelling is of the 120" casing described in previous reports. 

The quadrilateral shell element was used to allow modelling with a coarser 
mesh at a small penalty in solution accuracy. The geometric model used for 
the present analysis is shown in Figure 5.1 and the mesh is shown in Figure 
5.1A. 

LOADING 

The mechanical loading distribution used for this analysis' is shown 
on Figures 7.1 and 7.2. Note the radical difference between this pressure 
loading and the pressure loadings used in previous analyses. The total applied 
pressure is about 2.5 times less in this analysis than that used previously. 

The above mechanical loading was reacted by an interia loading as 
previously derived. For this analysis, the necessary constants have the following 
values : 
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P x = -9.2288 x 10 4 lb 
z c P x = -4.66293 x 10 7 lb-in 
V = 7.06297 x 10 4 in 3 
J = 3.51284 x 10 7 in 4 
I = 2.40335 x 10 10 in 5 

where 

P x = net externa] force in x direction 

Z c = cen ter of influence of external pressure 

V - Volume of shell 

= 1st moment of area of shell about y axis 
ly = 2nd moment of area of shell about y axis 

RESULTS 

a) Displacements 

The full load deformation mode (Figures 7.3 and 7.4) shows peak 
displacement occurring more toward the forward end as compared 
with the previous load distribution: at Z ^ 480" in the current 
analysis, and Z - 750 with the previous load case. Surprisingly, 
although the pressure spike is considerably more localized in 
this analysis, the deformation mode does not appear to be cor- 
respondingly localized. There is a possibility that the mesh 
used for this analysis was excessively coarse for such a loading: 
no convergence tests were performed for the new pressure dis- 
tribution, since it wa s felt that, based on previous experience 
with such elements, this mesh could capture sharp strain 
gradients. Nevertheless, convergence tests would be of value. 
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Figure 7.5 shows peak x-direction displacement on the load symmetry 
plane plotted as a function of full load fraction. The linearity of the re- 
sponse is clearly evident in this plot, the initial slope change being attributed 
to the constant internal pressure. Peak displacement was 3.73", compared to 
about 7.5" in the previous analysis at 45% of that prescribed full load. This 
difference is not surprising, since little localization was observed, and the 
net load in the present case is 2.5 times less than that in section 6. 

b) Stress 

Maximum stresses occur at the element integration point closest 
to the load spike and are as follows: 


Equivalent stress (psi) 

33000 

Hoop stress (psi) 

22000 

Axial stress (psi) 

16000 


The various stress components calculated are on the inside and outside 
of the casing contoured on Figures 7.6-7.11. Note especially in the two plots 
of equivalent J ^ stress (Figures 7.6 and 7.7) the similar stress pattern with 
the maximum stress occurring directly beneath the peak pressure, and with high 
stress gradients in that neighborhood. 

Compared with previous analyses stresses are overall much lower again 
because the pressure distribution used for this analysis is less severe than 
those used in previous analyses. 

CONCLUSIONS 

The present analysis suggests that there is no likelihood of geometric 
or material failure under the prescribed loading pattern and magnitude. This is 
in sharp contrast to the previously analyzed load cases (where the same finite 
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element model was used) where even with the 2 psi internal pressure, geometric 
failure occurred at 45% of the prescribed load magnitude. 

All analyses to date have been based on neglect of the stiffening 
effects of the casing segment joints. Certainly this may be regarded as a 
conservative assumption, so that a stiffened analysis of the current load case 
should not be necessary, unless load magnitudes are to be significantly in- 
creased. 
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SUMMARY OF RESULTS AND CONCLUSIONS 
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From the results of the earlier analysis it may be concluded that 
the doubly curved triangular shell elements are not suitable for analysis 
of very flexible shell structures with small thickness to diameter ratios such 
as the rocket motor casings analyzed here. ■ 

The quadrilateral doubly curved shell elements were found to be much 
better than the triangular elements and were flexible enough for analysis of 
the current problems. The non-linear results using the original tentative 
loadings on the structure predicted that the structure would fail in service 
with a 35% application of full load. It was found that an internal pressure 
of 2 psi provided an extra load bearing capacity of 10%. It would be useful 
to show that the effect of this internal pressurization varies linearly with 
pressure. 

The analysis with the more accurate load conditions show no signs 
of non-linearity up to full load. It is thus concluded that the structure is 
safe for service under these loads. The linear results obtained here were in 
agreement with results obtained with NASTRAN and BOSOR. This gives further 
confidence to the use of analysis for design of such shells subjected to splash 
down loadings. 
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APPENDIX 1: DESCRIPTION OF THE MARC-CDC PROGRAM 
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MARC-CDC is a non-linear general-purpose program which attempts to in- • 
corporate the latest in finite-element technology. The program consists of 
various modules, each of which represents a particular aspect of finite-element 
technology -- such as element geometry, material behavior, structural mechanics 
and matrix manipulation. 

The modules may be selected and used simultaneously with each other in 
a variety of combinations, resulting in a broad-based program which achieves 
the full potential of a general -purpose program. Aimed at moderate-size pro- 
blems, the program offers up to 10,000 degrees of freedom (depending on band- 
width size) and maintains a balance between efficiency and size. 

Designed to operate on CDC 6600 computer systems, MARC-CDC was de- 
veloped by the MARC Analysis Research Corporation of Providence, R.I. It is 
available in the United States through Control Data Corporation's CYBERNET 
services network. 

PROGRAM CAPABILITIES 

Because of its many capabilities, MARC-CDC is well suited for state- 
of-the-art analysis in applications such as pressure-vessel design, nuclear 
reactor design, and aerospace and automotive research. Among the many modules 
of MARC-CDC are the following libraries, programs, and sub-programs: 

Element Library -- a library of 25 elements, including curved beams 
with open and closed sections in one, two and three dimensiona; plates and 
shells; and solid two and three-dimensional i soparametric elements. These 
elements, when used in combination, have been selected so they can model any 
conceivable structure, depending upon size limitations of the program. 
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A special linkage feature enables the interconnection of shell and solid 
elements. All elements incorporated within the total MARC-CDC library have 
been constructed to permit the use of large displacement behavior and non- 
linear material behavior theories. 

Non-linear Behavior Library -- This library contains a selection of 
theories which cover the major material nonlinearities encountered in structural 
analysis. Creep analysis is effected by an initial strain procedure. Elastic- 
plastic analysis with either isotropic or kinematic hardening (in full nine- 
dimensional stress space) and temperature-dependent material properties is per- 
formed with an incremental (tangent modulus) procedure. A pressure-dependent 
yield stress is available for use with soil and other pressure-dependent materials 
(generalized Mohr-Coulomb behavior). 

Structural Theory -- Program control has been implemented to permit 
the combination of both the element library and behavior library in a series of 
piece-wise linear incremental analyses of every type of structure. Accordingly, 
analysis problems involving creep, buckling, combined nonl inear material and 
geometric behavior may be performed. 

Compatible Heat-Transfer Analysis -- A separate nonlinear finite-element 
heat transfer program generates compatible thermal data for use with the stress 
program. Typical problems include fusion (latent heat) effects, temperature- 
dependent thermal properties and nonlinear boundary conditions (such as radiation). 

Pre-Processor Programs -- MARCMESH-2D and MARCMESH-3D are available for 
the generation and bandwidth optimation of two and three-dimensional problems. 

With MARCSTRUCPL, input data may be plotted in planes or perspective. 

Restart and Automatic Load Incrementation -- These features have been 
developed during extensive experience with nonlinear analysis and result in an :• 
economy of computing time. The actual flow of an analysis can be seen by following 
the flow sequence given in Figure A.l. 
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ADDITIONAL FEATURES 

MARC-CDC also permits users to substitute their own subroutines to 
perform special tasks, thereby increasing the overall flexibility of the total 
program. These subroutines may cover a wide range of functions, from special 
input and output, to special creep laws and boundary conditions for thermal analysis. 
In addition, the MARC Corporation is currently developing many other modules and 
capabilities to further enhance the program. 

PROGRAM CAPABILITIES USED IN PROJECT 

CURVED TRIANGULAR SHELL ELEMENT - 
MARC-CDC Element 8 

Because we are concerned with the use of the curved triangular shell 
element - Element 8, we shall give a brief description of this element and 
its major features. This element is an isoparametric curved triangular shell 
element based on the Koiter-Sanders shell theory, which fulfills continuity 
requirements and represents rigid-body motions exactly. 

GEOMETRY 

The middle surface of the shell is defined by the equations: 
x = x(e ls e 2 ) 

y = y(e r e 2 ) 

where 

(x, y, z) are Cartesian coordinates. 

( 0 1 , ©g ) denote Gaussian coordinates on the middle surface e 1 t' ' c-'M 
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The domain of definition in the plane (o-j, e 2 ) is divided into a 
mesh of triangles which are mapped onto curved elements on the middle surface 
2. The actual middle surface is approximated by a smooth surface E which has 
the same coordinates (x-y-z) and the same tangent plane at each nodal point of 
the mesh. Practically, the mesh is defined by the Gaussian coordinates 
( e 1 ^ , 8«-) of the nodal points, and the surface E is defined by the values of 
According to the terminology of MARC-CDC, the coordinates are, therefore, the 
set: 


°2i 5 3x(p i )/3e 1 , ax(p i )/ae 2 

yfpj). sy(p i )/3Q 1> a y (p i )/ 30 2 
z(P i ). 9z(p i )/80 1 , 3z(p i )/39 2 

where 

x(p^) stands for x(e^., e 2 -j ) » (e-j ., 0 2 - ) being the coordinates of 
the node p . . 

DISPLACEMENTS 

There are nine degress of freedom for each nodal point p.. . These 
degrees of freedom are defined in terms of the Cartesian components of dis- 
placement u,v, and w, and rates of change with respect to the Gaussian 
coordinates. 
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u(p i ), autp.J/ae-,, au(p i )/ae 2 

“ i 

v(p 1 ) , av(p 1 .)/ae i# sv(p i )/3e 2 

V. 

—1 

w(p 7 - ) , awfp^/ae^ 3 w(p | J/ 39 2 
*i 

The displacements within an element are defined by interpolation 
functions. These interpolation functions g 2 ) are such that compati- 

bility of displacements and their first derivatives is insured between ad- 
jacent elements. Hence, for an element whose vertices are the nodal points 
p.{ , Pj and p^, the components u, v, w are defined as: 

u(e 1( s 2 ) .(u{, uT, uj) . <f(G r e 2 ) 
v ( e i , e 2 ) = (vl, vT, vj) . *(0 lf 0 2 ) 
w(e r 0 2 ) = (w{, wj, wj) . <,(o r 0 2 ) 

NUMERICAL INTEGRATION 

For this element, seven points of integration are used, with a rule 
which is exact for all polynomials up to the fifth order. 

LOADS 

Three different distributed loads are possible. The first case is a 
uniform load proportional to the surface area, positive when applied in the 
negative z-direction. The second case is a uniform normal pressure per unit 
of area on the surface. The third case is a non-uniform normal pressure per 
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unit surface area. In addition to uniform loads, concentrated loads may be 
applied at the nodes. 

LARGE DEFORMAT I ON ANALYSIS 

The large deformation analysis follows the Lagrangian description used 
in MARC-CDC [A1 ] . Note that only large deflection terms corresponding to the 
stretching strains have been introduced. This approximation is usually acceptable 
even for nonlinear buckling analysis. 

CURVED QUADRILATERAL THIN-SHELL ELEMENT - 
ELEMENT 4 

This is an isoparametric, doubly-curved thin shell element, using 
bi-cubic interpolation functions. The element is based on Koi ter-Sanders 
shell theory, fulfilling continuity requirements, and represents rigid body 
modes exactly when used as a rectangle in the mapped surface coordinate plane. 

The element contains no patching functions, so that it is restricted to qua- 
drilateral meshes with a maximum of four elements sharing one common node. 

However, the element is rapidly convergent in most problems which allow such 
a mesh. Note that any suitable surface coordinate systems may be chosen, so 
that the mesh need not be rectangular on the actual surface. 

GEOMETRY 

The element is isoparametric, so that the actual surface is interpolated 

1 2 

from nodal coordinates. The mesh is defined in the 0 -0 plane of surface 
coordinates. Then the actual surface is approximated by a surface defined by 
cubic interpolation on the interior of each, element based on the following 
set of 14 nodal coordinates: 


0 


1 



X, 


ax ax ay ay 
ao 1 ’ ae 2, y ’ ae'"* ae 2 * 


Z» 


2 2 2 
az az a x ay a z 

i 2 1 ? * l ? 1 ? 

30 ’ ao ’ 30 30 ’ 30 atr* 38 30 - 
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In most practical cases, the surface is definable as: 
x = x(e , © ) 

y = y(e\ e Z ) 

z = zfe 1 , e 2 ) 

Then the usual procedure is to define the mesh in the ©-e plane (as a rectangular mesh) 
by supplying the first two coordinates, o\ e 2 , at each node, through the COORDINATE 
input option. Then the remaining 12 coordinates are defined at each node through the 
use of the subroutine FXORD. 

The thickness of the element is defined in EGE0M1 . Note that when elements of 
differing thickness abut, tying must be used to avoid the imposition of improper 
continuity of membrane strain. 


DISPLACEMENT 


There are 12 degrees of freedom at each node: 

2 

3u 3u 3v 3v 3w_ 3w_ 3 u 

u ’ 30 1 ’ 30 2, v ’ 30 1 * 30 2, w ’ 30 1 5 30 2, 30 1 


these are 
3 2 v 


30 


2 1 2 

* 30 30 ’ 


3 2 W 

30 1 30 2 


where u, v, w are the Cartesian components of displacement. 

Displacement is interpolated by complete bi-cubic on the interior of an element, 
so that equality of the above nodal degrees of freedom at the coincident nodes 
of abutting elements ensures the necessary continuity required for thin shell 
theory. 


Note that fixed displacement boundary conditions should never be associated with 
all 12 degrees of freedom at each node, since 3 degrees of freedom must always 
determine middle surface (membrane) strains at the node. Care must therefore be 
exercised both in the specification of kinematic boundary conditions-- they must 
be fully, but not over-fully specified, and in the application of moments so that 
the generalized forces and the conjugate displacement multiply together to give 
a rate of mechanical work. 
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CONNECTIVITY SPECIFICATION AND NUMERICAL INTEGRATION 


The nodal point numbers of the element must be given in anti -clockwise 

12 12 
order on the 0 -0 plane, starting with the point (min 0 and 0 ). Thus, in 

Figure A-2, the connectivity must be given as i, j, k, 1. 



1 


k 



* « » 

*5.5 



i 


i 







FIG. A- 2 FORM OF ELEMENT 4 

The element is integrated numerically using 9 points (Gaussian quadrature). 

The first integration point is always closest to the first node of the element, 
then the integration points are numbered as shown in Figure A-2. Point 5 
(controid of the element in the 0-0 plane is used for stress output if 
NSTRESS=0 (or ALL POINTS option not flagged). 

LOADS 

Distributed Loads 

Three specifications of distributed loads are available with this element. 

These are described below. d below. 

Distributed Load Type 1 

This gives a uniform load, proportional to surface area, positive in the negative 
2 -direction (self-weight or snow load). 

Distributed Load Type 2 

Uniform normal pressure per unit area of surface. The magnitude of the pressure 
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is positive when applied in the negative _n direction where n_ = a_-| x is the 
normal to the surface: a^ is tangent to the positive 8-j line, a^ is tangent 

to the positive 02 line. 

Distributed Load Type 3 

This is a normal pressure per unit area of surface and is the same as type 2 except 
the magnitude is a function of position. The magnitude is defined in the user 
subroutine FORCEM: the required header cards are: 

SUBROUTINE FORCEM(PRESS ,TH1 ,TH2,NN,N) 

User Coding 

RETURN 

END 

with PRESS = magnitude of pressure at this point, defined by the user 

/ 

in this routine. 

\ 

TH1 Surface coordinates of the point, passed in for use in 
TH2 [ this routine. 

NN Integration point number. 

N Element number. 

This subroutine will be called once per integration point in all elements of type 
4 listed with load type 3. For such elements, the magnitude of load input in 
TRACTIONS, card series 5, will be ignored and the value defined in the user subroutine 
FORCEM will be used instead. Note that TH1 , TH2, NN and N must not be changed. 
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